4. Modulazione di ampiezza

Sommario

Esercizi sulla modulazione di ampiezza delle misure mediante analisi dello spettro e filtraggio offline.

1 Librerie

Quanto segue richiede questi pacchetti:

library(tidyverse)
library(plotly)
library(gsignal)

Se necessario installarli con install.packages() o mediante la GUI.

2 Funzioni di utilità

Per comodità definiamo alcune funzioni di utilità.

Generazione di un segnale sinusoidale con diverse armoniche, come definite in una tabella (pars) con le colonne w (ampiezza), f (frequenza) e phi (fase) per le varie frequenze (una riga per ogni armonica):

signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
  })
}

3 Simulazione del segnale

Costruiamo l’uscita di uno strumento di misura, la sua modulazione in ampiezza, il rumore interferente e grafichiamo gli spettri.

Una nota su plotly

Nel corso di statistica viene spiegato l’uso di GGPlot2. Quest’ultimo è studiato per l’efficacia e chiarezza di presentazione, per le quali è superiore a plotly.

È comunque possibile rendere i grafici GGPlot2 interattivi, passandoli alla funzione ggplotly().

Si noti tuttavia che avere molti grafici interattivi (e con molti punti) rallenta parecchio il caricamento della pagina nel browser: meglio evitarlo per i grafici “intermedi” e riservarlo ai grafici finali.

Definiamo alcuni parametri generici e la tabella delle armoniche:

# Parametri
fc <- 2000        # Frequenza di campionamento (Hz)
fm <- 10          # Frequenza del segnale di misura (Hz)
fr <- 15          # Frequenza del rumore (Hz)
fp <- 300         # Frequenza della portante (Hz)
m <- 0.8          # Indice di modulazione (0 < m < 1)
duration <- 10    # Durata del segnale (secondi)

# tabella delle armoniche del segnale di misura
pars_m <- tibble(
  w = c(1, 0.5, 0.3), 
  f = c(fm, 2*fm, 4*fm), 
  phi = c(0, 0, 0)
)

pars_m %>% knitr::kable()

w	f	phi
1.0	10	0
0.5	20	0
0.3	40	0

Generiamo il segnale di misura, perturbato da un rumore normale:

ym <- tibble(
  t = seq(0, duration, by = 1/fc),
  y = signal(t, pars_m),
  yn = y + rnorm(length(t), 0, pars_m$w[1]/10)
) %>% 
  mutate(
    f = 0:(length(t)-1)/max(t),
    fft = fft(y)
  ) 

ym %>%  
  mutate(intensity = Mod(fft) / length(t)*2,) %>% 
  ggplot(aes(x=f, y=intensity)) +
  geom_line()

Generiamo il rumore interferente con il suo spettro. Parametri delle armoniche corrispondenti al rumore interferente:

pars_r <- tibble(
  w = c(1, 0.5, 0.8), 
  f = c(fr, 3*fr, 9*fr), 
  phi = c(0, 40, 90)
)
pars_r %>% knitr::kable()

w	f	phi
1.0	15	0
0.5	45	40
0.8	135	90

Aggiungo il rumore interferente alla tabella ym:

ym <- ym %>% 
  mutate(
    y_in = signal(t, pars_r),
    fft_in = fft(y_in)
  )

Ottengo il grafico dello spettro:

ym %>% 
  select(f, fft, fft_in) %>% 
  pivot_longer(
    contains("fft"), 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ ifelse(.=="fft", "signal", "intf. noise")
  ) %>% 
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line()

Tidy data

Cerchiamo più possibile di utilizzare dati tidy, cioè una variabile per colonna, un’osservazione per riga. In questo modo realizzare grafici con serie multiple è particolarmente rapido (basta identificare la variabile di classificazione, cioè il nome della colonna che contiene la chiave di raggruppamento).

Per questo motivio, qui e di seguito useremo pivot_longer.

Aggiungo anche la modulazione:

# Aggiungo una mappa da nomi delle colonne -> nomi in legenda
fft_names=c(
  "fft" = "signal", 
  "fft_in" = "intf. noise", 
  "fft_m" = "modulation"
)

# Aggiungo la modulazione
ym <- ym %>% 
  mutate(
    y_m = (1 + m * yn) * cos(2 * pi * fp * t),
    fft_m = fft(y_m)
  )

# Grafico dello spettro
(ym %>% 
  select(f, contains("fft")) %>% 
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>%
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0, 500))) %>% 
  ggplotly()

Alternativa SENZA pivot_longer

Ridotto all’osso, l’ultimo grafico è (evito ggplotly per brevità):

ym %>% 
  select(f, contains("fft")) %>% 
  pivot_longer(-f) %>%
  mutate(
    intensity = Mod(value) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=name)) +
  geom_line()

Se non voglio usare pivot_longer() devo aggiungere tre colonne, una alla volta (stando attento a ripetere la formula correttamente: se scrivessi fft_m = Mod(fft_n) / length(t) * 2 otterrei un errore molto difficile da individuare!), con il modulo della FFT e nel plot devo aggiungere un layer alla volta (stando asttento a specificare correttamente il nome della serie come variabile color):

ym %>% 
  select(f, contains("fft")) %>% 
  mutate(
    fft = Mod(fft) / length(t) * 2,
    fft_in = Mod(fft_in) / length(t) * 2,
    fft_m = Mod(fft_m) / length(t) * 2
  ) %>% 
  ggplot(aes(x=f)) + 
  geom_line(aes(y=fft, color="fft")) +
  geom_line(aes(y=fft_in, color="fft_in")) +
  geom_line(aes(y=fft_m, color="fft_m"))

4 Effetto interferente sul segnale modulato

Adesso che abbiamo il segnale modulato separato in frequenza dal rumore possiamo ‘esporlo’ al rumore interferente.

Osserviamo prima la modulazione:

ym %>% 
  ggplot(aes(x=t, y=y_m)) + 
  geom_line()

Combiniamo segnale e modulazione:

ym <- ym %>% 
  mutate(
    y_m_in = y_m + y_in,
    fft_m_in = fft(y_m_in)
  )

Ora ym è così composta:

ym %>% head() %>% knitr::kable()

t	y	yn	f	fft	y_in	fft_in	y_m	fft_m	y_m_in	fft_m_in
0.0000	0.0000000	0.1418529	0.0	0.00e+00+0.0000000i	1.1213938	1.121394+0.0000000i	1.1134823	3.829540+0.000000i	2.2348761	4.950934+0.000000i
0.0005	0.1004060	0.1040175	0.1	2.10e-06-0.0132478i	1.1483848	1.121397-0.0073390i	0.6366972	4.564120+5.526609i	1.7850820	5.685517+5.519270i
0.0010	0.2000641	-0.0226319	0.2	8.30e-06-0.0265019i	1.0386500	1.121407-0.0146797i	-0.3034221	-3.960835-5.672208i	0.7352280	-2.839429-5.686887i
0.0015	0.2982363	0.3224292	0.3	1.87e-05-0.0397689i	0.8266720	1.121422-0.0220241i	-1.1963752	4.411363-12.118149i	-0.3697032	5.532785-12.140174i
0.0020	0.3942043	0.4672451	0.4	3.33e-05-0.0530552i	0.5637101	1.121445-0.0293739i	-1.1114244	-5.362537+4.257905i	-0.5477143	-4.241092+4.228531i
0.0025	0.4872785	0.5019278	0.5	5.21e-05-0.0663672i	0.3085893	1.121474-0.0367308i	0.0000000	4.785527+5.336739i	0.3085893	5.907001+5.300008i

ym %>% 
  select(t, y_m, y_m_in) %>% 
  pivot_longer(-t, names_to = "series", values_to = "value") %>% 
  ggplot(aes(x=t, y=value, color=series)) +
  geom_line() +
  facet_wrap(~series, nrow=2) +
  xlim(c(0, 0.5))

Osserviamo lo spettro:

fft_names["fft_m_in"] <- "mod. + intf."
ym %>%
  select(f, contains("fft")) %>% 
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>% 
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0,500))

Esercizio

Estrarre il segnale in uscita originale a partire da quello affetto da rumore bianco, modulato ed affetto da rumore interfente ovvero da ym$y_m_in

4.1 In un colpo solo

Nella parte precedente, la tabella ym è stata costruita gradualmente, allo scopo di illustrare la composizione del segnale simulato. Guardando le operazioni nel loro complesso, la procedura è in realtà più semplice:

tibble(
  t = seq(0, duration, by = 1/fc),
  y = signal(t, pars_m) + 
      rnorm(length(t), 0, pars_m$w[1]/10),    # segnale 
  fft = fft(y),                               # e sua FFT
  y_in = signal(t, pars_r),                   # rumore interferente
  fft_in = fft(y_in),                         # e sua FFT
  y_m = (1 + m * y) * cos(2 * pi * fp * t),   # modulante
  fft_m = fft(y_m),                           # e sua FFT
  y_m_in = y_m + y_in,                        # modulante + interferente
  fft_m_in = fft(y_m_in)                      # e sua FFT
) %>%
  # aggiungo la frequenza
  mutate(
    f = 0:(length(t)-1)/max(t),
  ) %>% 
  # prendo solo le colonne frequenza e fft
  select(f, contains("fft")) %>% 
  # riorganizzo la tabella
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>% 
  # aggiungo l'intensità
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0,500))

5 OMETTERE E FAR FARE COME ESERCIZIO

Passo 1: Filtraggio e Demodulazione

Estraiamo il segnale modulato mediante un filtro passa-alto. Di conseguenza è rimasto solo il segnale modulato che va quindi demodulato per ottenere una stima del segnale in uscita originario

cutoff <- 2 * max(pars_r$f) / fc  # Frequenza di taglio (2 volte la massima frequenza del rumore)
bf <- butter(4, cutoff, type = "high")
ym <- ym %>% 
  mutate(
    y_m_f = filtfilt(bf, y_m_in) # filtfilt non introduce ritardo
  )

Passo 2: Demodulazione

Moltiplicazione per la portante:

ym <- ym %>% 
  mutate(
    demod_raw = 2 * y_m_f * cos(2 * pi * fp * t) - 1
  )

Passo 3: Filtraggio Passa-Basso con filtro Butterworth

cutoff <- 2 * max(pars_m$f) / fc  # Frequenza di taglio (2 volte la massima frequenza del segnale informativo)
bf <- butter(4, cutoff, type = "low")
ym <- ym %>% 
  mutate(
    demod = filtfilt(bf, demod_raw)
  )

Confronto dei segnali:

(ym %>% 
  mutate(y_noisy = y + y_in) %>% 
  select(t, y, y_noisy, demod) %>% 
  pivot_longer(-t) %>% 
  ggplot(aes(x=t, y=value, color=name)) + 
  geom_line() +
  coord_cartesian(xlim=c(0, 0.5))) %>% 
  ggplotly()

Si noti come il segnale originario fosse affetto da rumore bianco, poi modulato e quindi sovrapposto ad un effetto interferente a tre armoniche. Alla fine del processo abbiamo ottenuto il segnale molto simile al segnale originario con una riduzione vistosa anche del rumore interferente